Geodynamic models of Indian continental flat slab subduction with implications for the topography of the Himalaya-Tibet region

The slab structure and high elevation of the Himalaya-Tibet region and their underlying mechanisms have been widely discussed. Many studies interpret a flat slab segment of Indian continental lithosphere located below the overriding plate, but interpretations of the northward extent of the flat slab differ substantially, with minimum estimates placing the boundary at the northern margin of the Himalaya (Indus-Yarlung Tsangpo suture), and maximum estimates placing it at the northern boundary of Tibet. In this study, we investigate for the first time if a flat slab segment of subducted buoyant Indian continental lithosphere below the Himalaya-Tibet region is geodynamically feasible and we quantify its northward extent, as well as its contribution to the high topography of the region. We conduct three large-scale fully-dynamic (buoyancy-driven) analogue experiments to simulate the subduction of the Indian continent. Our preferred, and geodynamically most feasible, model shows a continental flat slab extending northward up to ~ 320 km from the Himalayan thrust front, in agreement with recent estimates. Furthermore, it suggests that the positively buoyant flat slab segment of the Indian continent contributes some ~ 1.5–2 km to the high topography of the Himalaya-Southern Tibet region by providing an upward force to elevate the overriding Eurasian plate.


Geodynamic
The slab structure and high elevation of the Himalaya-Tibet region and their underlying mechanisms have been widely discussed.Many studies interpret a flat slab segment of Indian continental lithosphere located below the overriding plate, but interpretations of the northward extent of the flat slab differ substantially, with minimum estimates placing the boundary at the northern margin of the Himalaya (Indus-Yarlung Tsangpo suture), and maximum estimates placing it at the northern boundary of Tibet.In this study, we investigate for the first time if a flat slab segment of subducted buoyant Indian continental lithosphere below the Himalaya-Tibet region is geodynamically feasible and we quantify its northward extent, as well as its contribution to the high topography of the region.We conduct three large-scale fully-dynamic (buoyancy-driven) analogue experiments to simulate the subduction of the Indian continent.Our preferred, and geodynamically most feasible, model shows a continental flat slab extending northward up to ~ 320 km from the Himalayan thrust front, in agreement with recent estimates.Furthermore, it suggests that the positively buoyant flat slab segment of the Indian continent contributes some ~ 1.5-2 km to the high topography of the Himalaya-Southern Tibet region by providing an upward force to elevate the overriding Eurasian plate.
The Himalaya-Tibet mountain system in Asia, the highest mountain range in the world, has an average elevation of ~ 5 km (Fig. 1) 1 , formed mainly due to tectonic processes since the India-Eurasia continent-continent collision started at ~ 50 Ma [2][3][4][5][6][7][8][9] .The Himalayan mountain belt extends from the India-Eurasia collision zone plate boundary northward over a distance of ~ 200-320 km, until the southern margin of Tibet (Indus-Yarlung Tsangpo suture), while the Tibetan plateau extends northward over a distance of ~ 1000 km.Different conceptual models have been proposed to explain the formation of the high topography of the Himalaya-Tibet region, such as the underthrusting model 6,[10][11][12][13] , the crustal thickening model through localized 14,15 or distributed deformation 16,17 , and the lower crustal flow model 18 .Argand 13 proposed that the low-angle underthrusting of the Indian slab doubled the continental crustal thickness, which elevated the overriding plate above it.Others have also argued that most or all of the Himalaya-Tibet region has been underplated by Indian continental lithosphere using geological arguments 6 , or geophysical constraints [19][20][21][22] .However, other observational studies imply that the underthrusting of the Indian continent below the overriding Eurasian plate only reaches the northern margin of the Himalaya 12,[23][24][25] or a boundary located somewhere in Southern Tibet [26][27][28][29][30][31] .
Another point of scientific discussion concerns the possibility of continental subduction at the India-Eurasia convergent zone.Several researchers have proposed that subduction of the positively-buoyant Indian continent is not feasible unless there is significant eclogitization during subduction that increases the density of the continental crust [33][34][35] .Some modelling works 36,37 also considered that subduction of the entire crust of the Indian continent is infeasible and proposed that the upper crust has been scrapped off by the overriding plate during continental subduction leaving only the lower Indian crust and lithospheric mantle to subduct into the ambient mantle.However, the above studies did not consider the existence of laterally connected subduction zones, namely the Sunda subduction zone and Makran subduction zone, which also contribute to Indian continental subduction and northward motion, and East Asian deformation [38][39][40][41] .

OPEN
Department of Earth Sciences, Vrije Universiteit Amsterdam, Amsterdam, The Netherlands.* email: k.xue@vu.nlTherefore, geodynamic modelling is required to investigate if Indian continental subduction and flat slab formation are geodynamically feasible (i.e.dynamically feasible on Earth in terms of driving and resistive forces in the context of plate tectonics, subduction and mantle flow), when considering the laterally connected subduction zones (Sunda and Makran subduction zones), if and how the flat slab affects the topography of the Himalaya-Tibet region, and to provide quantitative insights on its horizontal extent.Earlier studies focusing on the subduction process of oceanic lithosphere have shown that subduction of buoyant features on the oceanic lithosphere (e.g.aseismic ridges, seamount chains, oceanic plateaus) can produce flat slab subduction below the base of the overriding plate [42][43][44][45] .Such a process can increase the surface topography of the overriding plate directly above the flat slab segment because of the higher average buoyancy of the oceanic lithosphere segment carrying the buoyant feature.Since the Indian continent is positively buoyant with a considerable size and seismic images also show underthrusting of the Indian slab segment at a shallow depth 46 , it is worthwhile to investigate whether subduction of the positively buoyant Indian continent contributes to the high surface topography of the Himalaya and Tibet.Our experimental models build on the upper mantle models of Bose et al. 38 , which demonstrated that lateral subduction zones (Sunda and Makran) are required for northward Indian continental motion and Indian indentation into the Eurasian plate, but which showed relatively limited continental subduction and no continental flat slab subduction.Our models will test the role of deep mantle subduction in driving continental subduction and the formation of a continental flat slab.
In this study, we present three 4-dimensional (3D space + time) buoyancy-driven (fully-dynamic) analogue experiments to investigate whether deep mantle subduction enhances continental subduction of the buoyant Indian lithosphere, if continental subduction can result in a flat slab directly below the base of the overriding plate lithosphere, and, if so, whether the flat slab contributes to the high surface topography above it.Since the amount of Indian continental subduction can be affected by the available slab pull and the resistance of the ambient mantle, we vary the mantle depth and lower-upper mantle viscosity ratio (η LM /η UM ) between the experiments: Exp. 1 only has an upper mantle reservoir with an infinitely high η LM /η UM at the 660-km discontinuity, Exp. 2 has a deep mantle reservoir with no viscosity step at the 660-km discontinuity (η LM /η UM = 1), and Exp. 3 has a deep mantle reservoir with an intermediate η LM /η UM (~ 8.6).

Methodology
We present three buoyancy-driven analogue models to simulate the India-Eurasia-Sunda collision-subduction zone in a big tank with a length of 180 cm and a width of 150 cm, which is filled with different amounts of glucose syrup, either homogeneous or layered, to represent the sub-lithospheric mantle reservoir (Fig. 2).Our models build on the experimental design presented in Bose et al. 38 , who presented buoyancy-driven, upper mantle, India-Eurasia collision-subduction experiments that either included or excluded lateral subduction zones.We improve on the experimental set-up of their most successful experiment (which includes the lateral Sunda and Makran subduction zones) by presenting experiments that include a lower mantle reservoir below the upper mantle.
A single layer of glucose syrup (density ρ UM = 1408 kg/m 3 and viscosity η UM = ~ 40 Pa s) with a thickness of 8.25 cm (scaling to 660 km in nature) is implemented in Exp. 1 to model the upper mantle reservoir, such that the rigid bottom boundary of the tank simulates the 660 km discontinuity with an infinitely high η LM /η UM .A thicker (18.125 cm, scaling to 1450 km in nature) glucose syrup reservoir is implemented in both Exps. 2 and 3.In Exp. 2, the glucose syrup is homogeneous with a density of 1408 kg/m 3 and a viscosity of ~ 40 Pa s to model a mantle reservoir with η LM /η UM = 1.In Exp. 3, the syrup reservoir is stratified with a density of 1396 kg/m 3 , viscosity of ~ 28 Pa s and thickness of 8.25 cm for the upper layer, representing the upper mantle, and a density of 1421 kg/m 3 , viscosity of ~ 240 Pa s and thickness of 9.875 cm for the lower layer, representing the lower mantle, giving a η LM /η UM of ~ 8.6.The value of η LM /η UM in Exp. 3 is slightly lower than the lower bound of the estimated ratio in nature (10-100) [47][48][49] .However, the resulting slightly lower resistance to subduction at the 660-km discontinuity is partly compensated for by the moderate lower-upper mantle density contrast (Δρ LM-UM = 25 kg/ m 3 ).We note that to dynamically scale our analogue experiments of subduction and collision, we choose to scale the density contrasts between the different layers and plates in our experiments (see below).When choosing this scaling approach, the absolute values of the densities are not important 50 .In our three experiments, these density contrasts are all the same, while the densities are slightly different.We used slightly lower densities for our Exp.3 compared to Exps. 1 and 2, because we had to reduce the viscosity of our sub-lithospheric upper mantle reservoir (syrup) in Exp. 3 by mixing it with water 51 .This was required to have the correct viscosity ratio between sub-lithospheric upper mantle and lower mantle.
For all three experiments, the overriding Eurasian plate and the oceanic lithosphere both have linear-viscous properties and are made of a mixture of high-viscosity, linear-viscous silicone and iron powder.The Indian continental crust consists of a mixture of 50 wt% of high-viscosity, linear viscous silicone, 39.62 wt% of high-viscosity, pink silicone rubber, and 10.38 wt% of low-viscosity, linear viscous silicone.The sizes of the plates are illustrated in Fig. 2 and a length scaling of 1 cm representing 80 km in nature is applied in this study.A wedge shape with an angle of ~ 30° is made at the leading edge of the Indian continent to simulate the passive margin.The overriding plate has a density of 1408 kg/m 3 for Exps.1-2 and 1396 kg/m 3 for Exp. 3, and is neutrally buoyant relative to the upper mantle material (ρ UM = 1408 kg/m 3 for Exps.1-2 and 1396 kg/m 3 for Exp.3).The oceanic lithosphere (ρ OL = 1508 kg/m 3 for Exps.1-2 and 1496 kg/m 3 for Exp. 3) has a density contrast of 100 kg/m 3 relative to the upper mantle material, which is slightly higher than that estimated in nature (80 kg/m 3 ) 52 to compensate for the ).A film of lubricant (made of 10wt% petrolatum and 90wt% paraffin oil) with a thickness of ~ 0.5-1.0mm is applied to the top of the subducting plate, in order to model a weak coupling at the subduction zone plate boundary interface, representing the hydrated uppermost crust and sediments at the subduction interface 53 .
Our models are dynamically scaled to nature according to the Stokes's settling law following earlier studies [53][54][55] ..In addition, the overriding Eurasian plate has a viscosity of ~ 2.8 × 10 4 Pa s (4.7 × 10 22 Pa s in nature), the subducting oceanic lithosphere has a viscosity of 3.0 × 10 4 Pa s (5 × 10 22 Pa s in nature), and the Indian crust has a viscosity of 2.5 × 10 4 Pa s (4.2 × 10 22 Pa s in nature).The scaling of time follows from the scaling of viscosity 50 , and in our experiments 1 s represents 8300 years in nature.We calculate the Reynolds number (Re) for our experiments following earlier studies 53 and the maximum values of Re are ~ 1.1 × 10 -3 , 6.9 × 10 -3 and 3.0 × 10 -3 for Exps.1-3, respectively.For all experiments, Re is much smaller than 1, which means our experiments have dynamic similarity with nature, reproducing a laminar flow regime with upstream-downstream symmetry 58 .
We use density contrasts between the plates and the mantle material to scale our models to nature.Accordingly, we need to apply a topographic correction factor (C Topo = ρ m_UM /ρ n_UM , where ρ m_UM and ρ n_UM represent the density of the upper mantle reservoir in the model and nature, respectively) to scale the topography in our models to that in nature, following Schellart & Strak 50 .Application of this method to our experiments indicates that 1 mm of topography in our models represents ~ 3.47 km in nature.
The experiments are initiated by pouring ~ 30 ml of upper mantle material on top of the first ~ 3 cm of the leading side of the subducting oceanic lithosphere, that forms an initial dip angle of ~ 30°.We calculate the topography of the overriding Eurasian plate with a stereo-photogrammetry technique with a subset size of 25 pixels and a step size of 6 pixels, following the approach as discussed in Chen et al. 59 .We record the slab geometry at the end of the experiments by taking photographs from the top after the overriding plate is manually pulled away to be able to view the geometry of the subducted slab from a top view perspective.The overriding plate was pulled away quickly to avoid the slab being deformed in a viscous manner, because at such high strain rates the response of the silicone mixture is mostly elastic 60 .

Results
Our experiments all start with rapid subduction of the leading oceanic lithosphere and the associated increase of the subducting plate velocity, trench velocity, convergence velocity and subduction velocity, and overriding plate extension.At the Sunda subduction zone, the trench motion is dominated by trench retreat from the beginning to the end of the experiments (Fig. 3).At the India-Eurasia collision zone, the motion of the plate boundary changes from retreat to advance and the convergence rate drops after reaching a maximum of ~ 14 cm/yr when the slab tip reaches the 660-km discontinuity for Exp. 1.In Exps. 2 and 3 the change from retreat to advance and the decrease in convergence rate after maxima of ~ 86 and ~ 26 cm/yr, respectively, occur when the Indian continent arrives at the plate boundary.In all experiments, the trench retreat at the Sunda subduction zone corresponds with a slab draping geometry, while the plate boundary advance at the India-Eurasia collision zone corresponds with a slab roll-over geometry in the experiments (Fig. 4).All experiments show continental subduction of the Indian continental lithosphere, but in Exp. 1 only minor continental subduction takes place (~ 1 cm scaling to ~ 80 km), while in Exps. 2 and 3 major continental subduction occurs (~ 12 cm scaling to ~ 960 km for Exp. 2 and ~ 7 cm scaling to ~ 560 km for Exp. 3) (Table 1, Fig. 4).

Topography
Before the onset of India-Eurasia collision, the frontal part of the overriding Eurasian plate is dragged down for all experiments (Fig. 3a,c,e, black and aqua lines in Fig. 5).This area is dragged down by a maximum of ~ 0.8, 0.4 and 0.3 mm (2.8, 1.4, 1 km) for Exps.1-3, respectively, from the beginning of the experiment (black lines in Fig. 5) to the onset of India-Eurasia collision (aqua lines in Fig. 5).The trench-normal extent of this area varies with time and per experiment, but generally ranges between ~ 10 and ~ 40 mm (scaling to 80-320 km).
From the onset of continent-continent collision and during the long-term Indian indentation, the frontal area of the overriding Eurasian plate is elevated in all experiments (Fig. 3b,d,f, Fig. 5), with a trench-normal extent of ~ 20, 60 and 35 mm (160, 480 and 280 km) for Exps.1-3, respectively.The least amount of uplift is observed for Exp. 1, which shows a maximum of ~ 0.5 mm (1.7 km) of uplift during the collision stage (i.e.maximum elevation difference between orange line and aqua line in Fig. 5a), and the elevation at the end of the experimental run is below the one observed at the start of the experiment by ~ 0.2 mm (compare orange and black lines).In contrast, Exps. 2 and 3 show significant overall uplift in the frontal area of the overriding plate (with respect to the zero elevation line), forming a bulge with a maximum height of ~ 2.1 mm (scaling to 7.3 km, orange line Fig. 5b) and 1.1 mm (scaling to 3.8 km, orange line Fig. 5c), and an average height of ~ 1.2 mm (scaling to 4.2 km) and 0.5 mm (scaling to 1.7 km), respectively.During the entire Indian indentation phase, the topographic bulge in Exps. 2 and 3 shows an overall movement towards the overriding plate, while the height of the bulge does not show any obvious change (compare blue and orange lines in Figs.5b,c).

Flat slab
After each experiment, the overriding plate is manually pulled away horizontally along the top surface to be able to view and investigate the geometry of the Indian continental slab.The geometry of the Indian continental slab is also visible and investigated from side view images (Fig. 4).For Exp. 1, only ~ 1 cm (80 km in nature) of the Indian continental lithosphere subducts below the overriding Eurasian plate and it dips steeply into the mantle (Fig. 4a).For Exps. 2 and 3, significantly more of the Indian continental lithosphere is subducted (12 cm, scaling to 960 km for Exp. 2, and 7 cm, scaling to 560 km for Exp. 3) (Table 1) (Fig. 4b,c).From this total length of subducted continental lithosphere, some ~ 44 (Exp.2) and 57% (Exp.3) closest to the trench underlies the overriding plate horizontally (Fig. 6), which we refer to as continental flat slab subduction, while the remaining ~ 56% (Exp. 2) and 43% (Exp.3) of the continental lithosphere dips steeply into the mantle.The trench-normal length of the continental flat slab area amounts to ~ 4.7-5.8cm (376-464 km in nature) and ~ 2.5-4 ± 0.5 cm (200-320 ± 40 km in nature) for Exps. 2 and 3, respectively.The western side of the flat slab area is somewhat narrower, compared to the eastern side for both experiments (Fig. 6).

High topography and continental flat slab subduction
Previous geodynamic models have shown that during trench advance due to oceanic subduction, the topography of the overriding plate near the trench is lower compared to the trailing part of the overriding plate and that this  Table 1.Measurements of the elevated topography (bulge) and flat slab length near the collisional plate boundary formed during long-term Indian indentation, and the amount of continental subduction.The length of the bulge and the flat slab are all measured along the mid axis of the Indian continent in the trench-normal direction.

Maximum height [mm] ([km]) Average height [mm] ([km])
Exp  ).An important observation is that the trench-normal extent of the continental flat slab and the topographic bulge are comparable (Table 1).On the other hand, in our Exp. 1 with only 1 cm (80 km) of continental subduction, both a topographic bulge at the surface and a flat slab are not observed.In addition, our models do not simulate crustal thickening due to orogenic wedge growth and the overriding Eurasian plate in our models is neutrally buoyant and so any thickening of this plate will not produce any topography.Therefore, we propose that the high topographic area of the overriding plate bordering the collisional plate boundary in our models results from the underlying continental flat slab.This can be explained by the positive buoyancy of the subducted Indian lithosphere, of which the subduction is driven by the slab pull from the Tethyan oceanic lithosphere in the upper and lower mantle, as well as the slab pull from the lateral Sunda and Makran subduction zones.The positive buoyancy of the Indian continental lithosphere (Indian crust + lithospheric mantle) promotes flat slab formation and provides an upward force to the base of the overlying Eurasian overriding plate, thereby forming a high surface topography.The trench-normal extent of the flat slab, the trench-normal extent of the bulge, and the height of the bulge are all higher in Exp. 2, compared to Exp. 3 (Table 1).The higher extent of the flat slab in Exp. 2 can be partly explained by the higher slab pull in Exp. 2 that resulted from the lower resistance of subduction at the 660-km discontinuity (low η LM /η UM ), compared to Exp. 3. The lower resistance at the 660-km discontinuity and large slab pull force promote faster subduction and cause more Indian continent to be dragged into the ambient mantle in Exp. 2 (~ 12 cm, scaling to 960 km), compared to Exp. 3 (~ 7 cm, scaling to 560 km).Considering that the continental lithosphere is positively buoyant, then more subduction of this lithosphere will promote a longer flat slab segment, and therefore the extent of the flat slab segment is also larger in Exp. 2 (5.3 cm long in the trenchnormal direction, scaling to 424 km), compared to Exp. 3 (4 cm long in the trench-normal direction, scaling to 320 km).The longer flat slab segment leads to a bigger bulge in the overriding plate above the flat slab in Exp. 2 (~ 32.5% longer than Exp.3).Additionally, the much faster subduction in Exp. 2 causes the mantle flow to be stronger, which results in a higher mantle wedge suction force, thereby promoting a longer flat slab segment and thus a larger topographic bulge compared to Exp. 3.
During the long-term Indian indentation and subduction, an obvious change of the bulge height is not observed in Exps. 2 and 3.This can be explained by the fact that the vertical forces applied to the overriding plate above the flat slab, which result from the buoyancy force of the Indian continent, do not vary much during continental subduction.The main change we observe is a northward displacement of the topographic bulge that can be attributed to the northward plate boundary advance.In our Exp.2, the slab pull is maximized, because there is no viscosity step at the 660 km discontinuity, producing a greater flat slab extent.However, in nature, the viscosity step is much higher than 1 and its estimated minimum is close to ~ 8.6 (as in our Exp.3).Thus, the extent of the flat slab, in the absence of other processes that could promote continental subduction and flat slab formation, should be lower than that estimated in our Exp.2, but more comparable to that observed in our Exp.3. Considering the above, our values of the maximum flat slab extent from Exp. 2 provide an upper bound as to what is geodynamically viable, even more so considering that the slab in Exp. 2 is continuous from the surface down to the lower mantle.Seismic tomography studies of the Himalayan region, however, imply that the slab is not continuous in many places along the collisional plate boundary 46,62,63 (Fig. 1a,b), likely due to slab detachment and break-off processes in the upper mantle.Such detachment and break-off will also lower the amount of slab pull 64,65 .
The importance of the conclusion of an upper bound for the extent of the continental flat slab as provided by Exp. 2 is that this upper bound is significantly lower compared to the values proposed by the above group 3. Furthermore, the results of our preferred experiment indicate that the continental flat slab can extend some 320 km northward from the MFT, suggesting that Indian underthrusting and continental flat slab subduction occur mainly below the Himalaya and part of Southern Tibet, in accordance with values reported in group 2 above.
In our experiments, the overriding plate experiences horizontal shortening due to northward indentation of India, with maximum shortening values of the order 80, 20 and 40% for Exps. 1, 2 and 3, respectively.The horizontal shortening and vertical thickening do not result in any topography, because our overriding plate is neutrally buoyant with respect to the sub-lithospheric upper mantle.Due to this neutral buoyancy, our experiments allow to directly quantify the contribution of the underthrusting of the Indian continental lithosphere to the uplift of the Himalaya-Tibet region.Our Exps. 2 and 3 present a bulge with an average elevation of ~ 4.2 and 1.7 km, which represent ~ 84% and 34% of the average elevation of the Himalaya (~ 5 km), respectively.Observational studies indicate a much higher crustal thickness reaching > 80 km beneath the Himalaya-Tibet orogen 66,67 compared to the average thickness of the continental crust in nature (~ 35 km 68,69 ), which suggests a considerable contribution of the crustal thickening to the high topography of the Himalaya-Tibet region.Since our Exp.3 has a viscosity step at the 660-km discontinuity that is the closest to nature, we suggest that the buoyant flat slab contributes as much as ~ 34% to the high topography of the Himalaya-Southern Tibet region, while other factors (e.g.crustal thickening) would contribute to the remaining part of the high topography.This contribution Vol:.( 1234567890

Comparison with previous models of flat slab subduction
To the best of our knowledge, these are the first reported buoyancy-driven geodynamic models of continental flat slab subduction.Flat slab subduction (with oceanic lithosphere) has been reported before in various modelling studies, both analogue and numerical, but these generally involved kinematic boundary conditions [43][44][45][70][71][72][73][74][75][76][77] . Only afew examples exist of buoyancy-driven modelling studies reporting (oceanic) flat slab subduction 78,79 .Some geodynamic modelling studies of flat slab formation at normal (oceanic) subduction zones have shown that the subduction of positively buoyant features (e.g.aseismic ridge, seamount chain, oceanic plateau) can lead to the formation of a flat slab [43][44][45]71 .Models from Espurt et al. 71 and Martinod et al. 45 both showed a flat slab and a high topography above the flat slab segment.However, the high topography (e.g. scalig to ~ 20 km in Espurt et al. 71 ) in their models is much higher compared to that in our models and in nature.This high topography could be caused by the high velocities, which resulted from external driving forces in their models.Our Exps.2-3 produce realistic topographies (maximum elevation reaching ~ 7.3 and 3.8 km, average elevation of ~ 4.2 and 1.7 km for Exps. 2 and 3, respectively) without an external driving force and this can be explained by the buoyancy-driven approach used in our models, that produce relatively lower velocities and forces during subduction and thus relatively lower topography.A recent buoyancy-driven analogue modelling study on the subduction of an aseismic ridge 80 did not produce a flat slab, but a lower slab dip angle at the ridge location was observed in one of their models, in which the ridge is the thickest (~ 20 km) and widest (200-400 km).The difference in outcome between their experiments (no flat slab) and our Exps. 2 and 3 (with flat slab) can be ascribed to various factors, including a much bigger size of the buoyant feature (Indian crust, scaling to 48 km thick and 2000 km wide), a much higher slab pull force because of the much deeper subduction of the slab into the lower mantle, and the existence of a much larger lateral subduction segment (Sunda) in our experiments.Earlier geodynamic models of the India-Eurasia collision have been presented, both analogue 81,82 and numerical 83 , but these have not reproduced continental flat slab subduction.Replumaz et al. 81 and Pitard et al. 82 presented India-Eurasia collision-continental subduction experiments that involved an external velocity boundary condition (piston) and which excluded lateral subduction zones.The absence of continental flat slab subduction in these experiments might be due to the relatively narrow width of the convergent plate boundary, scaling to only 1200 km in both studies (much narrower than in the current study, where it scales to 6400 km). Such a nrrow width would cause mantle wedge suction forces to be relatively small 78,84 , thereby preventing a flat slab to form.Pusok & Kaus 83 presented models with and without an external driving force for Indian continental subduction. A igh topography in front of the collisional boundary was only observed in their models with an external driving force, which was due to the crustal thickening, but not in their fully-dynamic model, which was driven only by the negative buoyancy of the oceanic slab. Ths could result from the insufficient subduction of Indian continent in their fully-dynamic models, thereby inhibiting a flat slab to form.This insufficient amount of Indian continental subduction could be explained by two main differences between their and our study.The first is the relatively lower negative buoyancy force in front of the Indian continent, which was caused by both the smaller mantle thickness implemented in their model (1000 km in their work vs. 1450 km in our models) and shorter Tethyan subducted segment before Indian continental subduction (< 1000 km in their work vs. ~ 1200 km in our models).The other difference is that their model implemented a narrower (~ 1500 km measured along the trench) Sunda subduction zone compared to that in nature and in our models (~ 3600 km), which diminished the role of the lateral oceanic subduction zone in driving continental subduction.

Conclusion
Our self-consistent, buoyancy-driven geodynamic experiments have demonstrated that large-scale underthrusting/continental flat slab subduction is physically viable in a continental subduction-indentation-subduction setting such as found at the India-Eurasia-Sunda convergent zone.Furthermore, our models provide new insights and understanding of the slab geometry at the India-Eurasia continental subduction zone and the high topography of the Himalaya-Tibet region.Our models demonstrate that a lower η LM /η UM promotes the subduction of positively buoyant continental lithosphere, while an increased amount of positively buoyant continental lithospheric subduction facilitates the formation of a flat slab.The positively buoyant flat slab that forms below the overriding plate provides an upward force to the overlying plate, forming a high surface topography at the leading edge of the overriding plate overlying the flat slab.Our models suggest that part of the high topography of the Himalaya and southernmost Tibet, on average possibly some 1.5-2 km, is related to the underthrusting of the positively buoyant Indian continental lithosphere directly beneath the base of the overriding Eurasian plate.Our self-consistent models provide a geodynamic constraint on the upper limit of the northern extent of Indian continental flat slab subduction, which is of the order 320 + /-40 km from the Main Frontal Thrust.

Figure 1 .
Figure 1.(a,b) Seismic tomography across the India-Eurasia collision zone from the global P-wave seismic tomography model of Amaru 32 .(c) Topographic map of the India-Himalaya-Tibet region illustrating the location of the tomography profiles (P1 and P2) shown in (a) and (b) and the topography profiles (AA' and BB') shown in (d) and (e).(d,e) topographic profiles across the India-Eurasia collision zone.Green arrows indicate the approximate downdip limit of the Indian continental flat slab as implied by the seismic tomography sections in panels (a) and (b).

Figure 2 .
Figure 2. Schematic model setup.Experiment 1 involves only upper mantle material with the rigid bottom of the tank simulating an infinite viscosity step at the 660-km discontinuity.Experiment 2 involves a homogeneous mantle layer with a thickness that represents 1450 km in nature and excludes a viscosity step at the 660-km discontinuity (η LM /η UM = 1).Experiment 3 incorporates an upper mantle layer with a thickness of 8.25 cm (660 km) and lower mantle layer with a thickness of 9.875 cm (790 km), with a viscosity step at the 660-km discontinuity (η LM /η UM = ~ 8.6).Two cameras using a stereo-photogrammetry technique are used to compute the overriding Eurasian plate topography from the top.
1 and the homogeneous upper and lower mantle in Exp. 2 have a viscosity of 40 Pa s, scaling to 6.73 × 10 19 Pa s in nature.Exp. 3 has a viscosity of ~ 28 Pa s (4.7 × 10 19 Pa s in nature) in the upper mantle and ~ 240 Pa s (4 × 10 20 in nature) in the lower mantle.The scaled viscosity values are all within the estimated range of 10 19 -10 21 Pa s for nature https://doi.org/10.1038/s41598-024-52709-wwww.nature.com/scientificreports/

Figure 3 .
Figure 3. Topographic maps of a sub-area of the overriding plate located at and near the India-Eurasia collisional plate boundary.The location of the sub-area is shown in the inset in the right of each panel.The white arrows represent the location of the topographic profiles in Fig. 5.The dotted lines in panels d and f indicate the down-dip limit of the continental flat slab.

Figure 4 .
Figure 4. Photographs of the three subduction-collision-continental subduction experiments in a late stage of the experimental evolution for (a) Exp. 1 (upper mantle experiment, (b) Exp. 2 (deep homogeneous mantle experiment), and (c) Exp. 3 (layered mantle experiment).The photographs show the mantle reservoir, the slab structure (oceanic slab is dark grey-black, Indian continental crust is light pink with thin black passive grid lines) and the bottom of the overriding plate.All experiments show a slab rollback and slab draping structure for the Sunda slab and a slab rollover structure for the Tethyan + Indian slab.Exps. 2 and 3 show Indian continental flat slab subduction.Photographs in (a) and (c) by K. Xue.Photograph in (b) by W.P. Schellart.
www.nature.com/scientificreports/low topography mainly results from the vertical component of the shear force at the subduction interface61 .Contrastingly, subduction of the buoyant Indian continental lithosphere during collisional plate boundary advance leads to a bulge (Fig.7) with high topography next to the plate boundary in our Exp. 2 (maximum height scales to ~ 7.3 km and average height scales to ~ 4.2 km) and Exp. 3 (maximum height scales to 3.8 km and average height scales to ~ 1.7 km

Figure 5 .
Figure 5. Diagrams showing the topographic profiles of the overriding Eurasian plate north of India at four times for Exps.1-3.(a-c) Start of subduction (start of the experiment) (black lines), start of continental collision and indentation (aqua lines), advanced stage of collision and indentation (blue lines), and final stage of the collision and indentation (end of the experiment) (orange lines) for the three experiments.The position of the topographic profiles is indicated by white arrows in Fig. 3.

Figure 6 .
Figure 6.Photographs of the continental flat slab for (a) Exp. 2, and (b) Exp. 3 which were taken after the overriding plate was manually pulled away to be able to view the subducted slab from above.Photograph in (a) shows a top-view perspective, while photograph in (b) shows a 3D perspective view.The solid line indicates the collisional plate boundary, while the dashed line delineates the downdip limit of the flat slab, which has a (sub)horizontal orientation and lies directly below the overriding plate, but not below the top surface of the upper mantle reservoir.(c) Photograph with a 3D perspective view of Exp. 2 at the end of the experimental run showing the elevated topography in the overriding plate at the India-Eurasia collision zone plate boundary.(d) Line drawing of photograph in (c).The white arrows in (d) indicate the direction of the plate boundary motion.Photographs in (a) and (c) by W.P. Schellart.Photograph in (b) by K. Xue.

Figure 7 .
Figure 7. Schematic cross-section illustrating the underthrusting of the experimental Indian continent below the experimental overriding Eurasian plate in our Exps. 2 and 3, forming a topographic bulge above the flat slab segment.The length scale of the bulge and flat slab refers to Exp. 2.
13 the high topography in Exp. 3 is limited to a region within ~ 320 km and located north of the collisional plate boundary.As such, it does not reproduce, nor does it provide, an explanation for the high elevation of the Tibetan plateau northward of this ~ 320 km boundary.Nevertheless, our models do indicate that underthrusting and continental flat slab subduction below the entire Tibetan Plateau are dynamically not viable, and thus our models do not support the conceptual models from e.g.Argand13and Powell & Conaghan 6 in which the high topography of the Tibetan Plateau results entirely from this underthrusting and continental flat slab subduction.